#include <cstdlib>
#include <stdlib.h>
#include <time.h>
#include <stdio.h>
#include <iostream>
#include <math.h>
#include <string.h>
#include <matrix.h>


#include "EllipticIntegrals.h"
#include "Coils.h"
#include "Psi.h"
#include "GreenFunction.h"
#include "Grid.h"
#include "Exception.h"
#include "Log.h"
#include "Constants.h"
#include "DynamicArray.h"
#include "Parameters.h"
#include "Jphi.h"
#include "Error.h"
#include "CombAlg.h"

#ifndef _NO_NAMESPACE
using namespace std;
using namespace math;
#define STD std
#else
#define STD
#endif

#ifndef _NO_TEMPLATE
typedef matrix<float> Matrix;
#else
typedef matrix Matrix;
#endif

#ifndef _NO_EXCEPTION
#define TRYBEGIN()	try {
#define CATCHERROR()	} catch (const STD::exception& e) { \
						cerr << "Error: " << e.what() << endl; }
#else
#define TRYBEGIN()
#define CATCHERROR()
#endif

Log* plog;
Parameters* ppmts;




Psi* psi;

float** psipl_1;
float** psipl_m;
float** psipl_tmp;

float** psist;
float** psitotal_1;
float** psitotal_m;
float** psitotal_tmp;
float** psipl_on_g; // for managing psi values on plasma boundary


float** psipl_s;
float** psipl_rs;
float** psipl_on_g_s;
float** psitotal_s;




float** jphi_0;
float** jphi_m;
float** jphi_tmp;

float** jphi_s;

Error* err;

int main(int argc, char** argv) {


    plog = new Log(LOG_FILE_NAME, PATH);
    plog->writeMsg("--- start program ---");

    ////////////////////////////// parameters ///////////////////////////////////////
    plog->writeMsg("initializing parameters");
    ppmts = new Parameters();


    ppmts->setCoilDataFile("CoilData.txt");
    ppmts->setLeftR(0.05);
    ppmts->setRightR(0.25);
    ppmts->setUpZ(0.2);
    ppmts->setDownZ(-0.2);
    ppmts->setNumberR(50);
    ppmts->setNumberZ(50);

    ppmts->setAlpha(0.5);
    ppmts->setBeta(0.73);
    ppmts->setBetaPl(0.2);
    ppmts->setPsiG(0.0005);
    ppmts->setI(50000);
    ppmts->setPrecision(0.0001);

    ppmts->setPlasmaLeftR(0.1);
    ppmts->setPlasmaRightR(0.2);
    ppmts->setPlasmaUpZ(0.1);
    ppmts->setPlasmaDownZ(-0.1);

    plog->writeMsg("initializing parameters finished");
  
    plog->writeMsg("reading coils data");
    Coils Gutta(ppmts->getCoilDataFile(), 4);
    plog->writeMsg("reading coils data finished");

    plog->writeMsg("initializing grid");
    Grid omega(ppmts->getLeftR(),ppmts->getRightR(), ppmts->getUpZ(), ppmts->getDownZ(), ppmts->getNumberR(), ppmts->getNumberZ());
    plog->writeMsg("initializing grid finished");


    plog->writeMsg("initializing psist");
   
    psi = new Psi(&omega, &Gutta);
    psist = psi->st_grid();
    plog->writeMsg("initializing psist finished");

    ////////////////////////////// jphi_0 ///////////////////////////////////////
    plog->writeMsg("initializing jphi_0");
    Jphi* jphi;
    jphi = new Jphi(&omega);
    jphi_0 = jphi->on_grid_initial();
    plog->writeMsg("initializing jphi_0 finished");

    ////////////////////////////// psipl_1 ///////////////////////////////////////
    plog->writeMsg("initializing psipl_1");
    psipl_1 = psi->pl_grid(jphi_0, psipl_1);

    plog->writeMsg("initializing psipl_1 finished");

    ////////////////////////////// psitotal_1 ///////////////////////////////////////
    plog->writeMsg("initializing psitotal_1");
    psitotal_1 = psi->total_grid(psipl_1, psist, psitotal_1);

    psipl_on_g = psi->get_psipl_on_g(jphi_0);
    plog->writeMsg("initializing psitotal_1 finished");

    ////////////////////////////// error ///////////////////////////////////////
    plog->writeMsg("initializing Error object");

    err = new Error(&omega);
    plog->writeMsg("initializing Error object finished");
    //////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////

    plog->writeMsg("beginning direct caclculations");
    float e = ppmts->getPrecision();


    float miscal = 0;
 
    int num_r = omega.getNumberR();
    int num_z = omega.getNumberZ();

    CopyDynamicArray(jphi_0, jphi_m, num_r, num_z);
    CopyDynamicArray(psipl_1, psipl_m, num_r, num_z);
    CopyDynamicArray(psitotal_1, psitotal_m, num_r, num_z);



   
   
    int k = 0;
    plog->writeInFile(JPHI_FILE_NAME, 0, jphi_m, &omega);
    plog->writeInFile(PSIPL_FILE_NAME, 0, psipl_m, &omega);
    plog->writeInFile(PSIPL_ON_G_FILE_NAME, 0, psipl_on_g, &omega);
    plog->writeInFile(PSITOTAL_FILE_NAME, 0, psitotal_m, &omega);


    if (true)
    {
      do
      {
        k+= 1;
        cout << "\n k: " << k << " miscal: " << miscal;
        plog->writeInFile(JPHI_FILE_NAME, k, jphi_m, &omega);
        plog->writeInFile(PSIPL_FILE_NAME, k, psipl_m, &omega);
        plog->writeInFile(PSIPL_ON_G_FILE_NAME, k, psipl_on_g, &omega);
        plog->writeInFile(PSITOTAL_FILE_NAME, k, psitotal_m, &omega);
        CopyDynamicArray(psipl_m, psipl_tmp, num_r, num_z);
        CopyDynamicArray(jphi_m, jphi_tmp, num_r, num_z);

        jphi_m = jphi->on_grid(psitotal_m, jphi_m);
    
        psipl_m = psi->pl_grid(jphi_m, psipl_m);
        psitotal_m = psi->total_grid(psipl_m, psist, psitotal_m);
        psipl_on_g = psi->get_psipl_on_g(jphi_m);
        miscal = err->getValue(psipl_tmp, psipl_m);
        
    } while (miscal > e && k < 500);
    
    plog->writeInFile("psitotal_end", psitotal_m, &omega);
    plog->writeInFile("psist", psist, &omega);
    psi->bound(psitotal_m);
    psipl_on_g = psi->get_psipl_on_g(jphi_m);
      plog->writeInFile("psipl_on_g_end", psipl_on_g, &omega);
    }


    plog->writeMsg("finished direct caclculations");

    ////////////////////////////////////////////////////////////////////////////
    // GRID COMBINED METHOD
    ////////////////////////////////////////////////////////////////////////////
    plog->writeMsg("beginning combined caclculations");
    ReadDynamicArray("res/jphi_end.txt", jphi_s, num_r, num_z);
    ReadDynamicArray("res/psipl_end.txt", psipl_s, num_r, num_z);

    CopyDynamicArray(psipl_on_g, psipl_on_g_s, num_r, num_z);

    plog->writeInFile("jphi_s", jphi_s, &omega);
    plog->writeInFile("psist", psist, &omega);
    float h_r, h_z;
    float ri, zj = 0;
    printf("\nplasma: %d %d %d %d", psi->left_r_g, psi->right_r_g, psi->up_z_g, psi->down_z_g);
    h_r = psi->h_r;
    h_z = psi->h_z;
    int rows, columns, n, m, dim = 0;
    m = psi->m;
    n = psi->n;
    CombAlg* calg2;

    calg2 = new CombAlg(&omega, psist, n+1, m+1);
    
    Matrix jphi_0 = calg2->jphi_s;
    Matrix psipl_g_0 = calg2->get_psipl_on_g(jphi_0);
    Matrix psipl_l = calg2->solve_sor2(psipl_g_0, jphi_0);
    Matrix jphi_1 = calg2->get_jphi(psipl_l);
    Matrix psipl_g;
    Matrix psipl_temp_l;
    psipl_temp_l.SetSize(n+1,m+1);
    int s = 0;
    //main loop
    if (false) {
        do {
            s += 1;
            psipl_l = calg2->solve_sor2(psipl_g_0, jphi_0);
            jphi_0 = calg2->get_jphi(psipl_l);
            psipl_g_0 = calg2->get_psipl_on_g(jphi_0);
            miscal = err->getValue(psipl_temp_l, psipl_l);
            cout << "\n s: " << s << " miscal: " << miscal;
            psipl_temp_l = psipl_l;

            plog->writeInFile(JPHI_FILE_NAME, s, jphi_0);  
            plog->writeInFile(PSIPL_ON_G_FILE_NAME, s, psipl_g_0);
            plog->writeInFile(PSIPL_FILE_NAME, s, psipl_l);
        }
        while (miscal > e*0.1 && s < 500);
        }
   
    delete(jphi);
    delete(psi);
    delete(ppmts);
    delete(err);
    delete(plog);
    delete(calg2);
    return 0;
}

